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Abstract 

C We introduce a technique of time series analysis, potential forecasting, which is based on 

I dynamical propagation of the probability density of time series. We employ polynomial 

+S coefficients of the orthogonal approximation of the empirical probability distribution 

^ and extrapolate them in order to forecast the future probability distribution of data. 

^ The method is tested on artificial data, used for hindcasting observed chmate data, and 

then apphed to forecast Arctic sea-ice time series. The proposed methodology com- 
^ pletes a framework for 'potential analysis' of climatic tipping points which altogether 

^ serves anticipating, detecting and forecasting climate transitions and bifurcations using 

, ^ several independent techniques of time series analysis. 

^ Keywords: potential forecasting, potential analysis, time series analysis 

> 

O 
On 
O 

Many dynamical systems in general, and geophysical subsystems in particular, lack 
^ analytical deterministic descriptions with fully developed physical models, being rep- 

^ resented mainly by recorded time series. At the same time such systems may be of 

^ great public interest and societal impact, such as the current climate change with rising 

• ^ temperature records around the globe. 

^ In these circumstances, powerful research tools may be provided by time series 

d analysis based on simple generalised stochastic models, where uncertain or unknown 

variables can be represented by stochastic components. 

In previous papers [HI 0, E], we have developed several time series techniques for 
anticipating and detecting tipping points in trajectories of dynamical systems, with 
applications in climatology. The modified degenerate fingerprinting [6] was introduced 
for early warning of critical behaviour in time series to allow one to anticipate an 
upcoming bifurcation or transition (climate tipping points ^). In order to distinguish 
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between possible transition and bifurcation we developed potential analysis [11 E], which 
is a method for deriving the number of system states under assumption of quasi- 
stationarity of data subset. If both techniques give indication of dynamical change, 
this denotes a genuine bifurcation. If modified degenerate fingerprinting indicates a 
change but potential analysis does not, this may mean a transition rather than a 
bifurcation, with no changes in the underlying system potential. 

In this paper, we develop the methodology further, so that the we become able 
not only anticipate and detect, but also forecast the time series dynamics. The skill 
of such a forecast will depend on several factors, in particular, whether the upcoming 
change will be gradual or abrupt, at what rate it will be happening and how the scaling 
properties of the stochastic component may change with time. 

Here we outline the methodology, test it on artificial data, in several hindcast case 
studies, and provide a forecast of the dynamics of Arctic sea-ice extent in the nearest 
future. 

2. Methodology 

2.1. Potential analysis as the basis of the method 

We consider a simple stochastic model with a polynomial potential U as an approx- 
imation of the system dynamics 



where z is the time derivative of the system variable z{t) (time series of an observed 
climatic variable), t] is Gaussian white noise of unit variance and a is the noise level. 
In the case of a double-well potential, it is a polynomial of 4th order: 



According to the Fokker-Planck equation for the dynamic evolution of the proba- 
bility density function p{z,t), 




(1) 






its stationary solution is given by [T] 




(3) 



The potential can be reconstructed from time series data of the system as 




(4) 
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which means that the empirical probabihty density pd has the number of modes cor- 
responding to the number of wells of the potential. 

This simple approximative approach allowed us to reconstruct the system potential 
of various climatic records (see [8]). It works with remarkable accuracy for data subsets 
of length as short as 400 to 500 data points, demonstrating above 90% rate of accurate 
detection, as was shown in an experiment with artificial data. For data subsets of 
length above 1000 points it correctly detects the structure of the potential with rate 
above 98% [8]. 

Here we develop the potential method beyond its detection capability, such that we 
are able to forecast the behaviour of a time series on the basis of its potential. To that 
effect, we introduce an extrapolation technique that would use the potential structure 
of the time series with linear extrapolation of the coefficients of the approximating 
polynomials. To reduce the biases introduced during various stages of the potential 
analysis (due to kernel distribution approximation, further logarithmic transformation, 
noise estimation, and finally polynomial fits), we use the empirical probability density 



rather than its logarithmic transformation, the potential (see Eq. |2.1 ). Moreover, 
unlike in previous work [8j, we use not the 2N potential coefficients (where N is the 
number of potential wells) but the coefficients of the approximation of the empirical 
probability density by a finite Chebyshev polynomial series (following the approach 
of [12j). Chebyshev approximation has an advantage of being near optimal, and already 
lOth-degree approximation in most cases of observed time series provides an accurate 
fit with low values of error function. 

2.2. Approximation of the probability density 

To approximate the empirical probability density, we use the orthogonal (in the 
interval [—1, 1]) Chebyshev polynomials of the first kind: 

To(x) = 1, 
Ti(x) = X, 

(5) 



T„+i(x) = 2xTn{x) - r„_i(x). 

The polynomial Tn{x) has n zeros in the interval [—1, 1] at points 

/7r(A;-l/2)\ 

X = cos , k = 1,2, . . . , n. 

\ n J 

The approximation of a function f{x) can be done by using a truncated (finite N) 
sum of the following form 

fi^)= I yckTkix) 1 -^co, (6) 



^CfcTfc(x) j - -Co, 

\k=0 J 
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where c„ are the coefficients obtained by discrete cosine transform of the vector of 
nonuniformly spaced samples of the considered function over the samphng grid 

X = cos ( — ^—^ 1 , k = 1,2 n. 



n 

For an arbitrary interval [a, h] it is necessary to transform variables as 

_ X- 1/2(6 + a) 
^~ 1/2(6-2) ■ 

A good example of the above calculations is given in |T2]. 

When decomposition (|6]) is obtained according to the particular time series prob- 
lem to be analysed, the resulting polynomial is expanded, thus producing the final 
coefficients 

f{^)=[Y.^kx\ (7) 

\A:=0 / 



2.3. Linear extrapolation of the coefficients and forecast time series 

We consider the approximation of the empirical probability density using Cheby- 
shev polynomials Tq, ... ,Tio. When the approximating polynomial is derived, the 
decomposition coefficients Ck (Eq. 7) are linearly extrapolated using a set of preceding 
values. The interval of these, as well as the extension of the extrapolated interval can 
be varied according to the particular time series to be analysed. 

Once a new probability density is calculated, we generate a forecast time series using 
rejection sampling algorithm (see, for instance, [2]). This provides an artificial series 
with the prescribed distribution, but this may be not enough for obtaining a realistic 
forecast time series, because the ordering of the series (and hence scaling properties 
like long-term memory) should be reconstructed according to the initial data. For 
this purpose, we apply so-called "sorting" of time series, i.e. arranging its values in 
the same order as in the initial data (before the forecast started), thus reproducing 
realistic correlations (because their distributions are already very similar due to the 
extrapolation of probability density). Sorting is a simple numerical algorithm which 
uses ranking of the values of two series, initial subsample and forecast subsample. 
However, this should be done with care, especially in data with seasonality: if there is 
a seasonal trend, it is very important to sort the forecast series according to observed 
data at the same date of the year, so that further seasonal variability would be adequate. 
This is achieved by going back along the series with step equal to the seasonality period 
(365 for daily data or 12 for monthly data). Since certain years may be anomalous in 
fluctuations (due to internal variability in the system), the initial data used for sorting 
may be an average over several years starting from the same date in a year (for instance, 
March 1st in several consequtive years). This average is then used to sort the forecast 
series starting on March 1st and projecting into future. 
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2.4- Uncertainties and applicability; criteria of performance 

In many cases of abrupt of highly nonlinear dynamics the linear extrapolation of 
the decomposition coefficients may produce the empirical distribution with large devia- 
tions, especially in case of non-stationarity of the data. We attempted bootstrapping of 
the decomposition of coefficients according to [H] . Based on bootstrapping techniques, 
it is possible to consider blocks of data in a chosen subset x of size 

^/{Q)ai{x -x) \ 
1 — (ai(x — x))^ / ' 

where W is the window length, ai is the lag-1 autocorrelation, x is the mean value of the 
subset X. In the case of nonstationary data, when the probability distribution varies 
within the data subset, bootstrapping provides estimates of the partial probability 
distibutions, which may deviate from the average quite significantly. In practical terms, 
applying bootstraping for estimation of the decomposition coefficients in non-stationary 
data provided worse results in the considered samples, with the forecast time series of 
poorer skill than the single-estimated probability density functions. A possible solution 
to this could be modification of the bootstraping algorithm, where instead of mean value 
removal a more sophisticated detrending is applied. We plan to adapt block bootstrap 
methods [H] (Chapter 3 therein) for that means in a future study. 

Furthermore, the important parameter that affects the skill of the forecast is the 
extrapolation period. The skill of the forecast drops with increase of its value. 

Obviously, in the case of abrupt changes using linear extrapolation of coefficients 
may prove unsatisfactory, and our method in those cases may be not applicable. The 
best results are obtained when the data undergoes gradual dynamic change and the 
forecast horizon is within 100 time units (which means 3 months for daily data and up 
to 8 years for monthly data). 

To assess the skill of the forecast, we used several techniques widely applied in 
modelling community for comparison of observed and modelled data, for instance, in 
hydrology 0: 




Daily Root Mean Square 



Nash-Sutcliffe efficiency 



Percent bias 



N 



X 100, 



where Qo and are the observed and modelled fluxes, respectively, i = 1, . . . , N is 
time index (daily flux), q is the mean value of the series q. It is easy to see that for 
two identical time series DRMS=0, NS=1, and %bias=0, and any deviation from those 
values would indicate the difference between the modelled and observed time series 
pointwise. 
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3. Tests of artificial and cHmate data 



3.1. Artificial data 

We considered several simulated time series: an artifical data where the potential 
is varying from single- well to double- well and back several times (Fig. [T]), double- well- 
potential data with decreasing noise level (Fig. |2]) and artificial data bifurcating from 
one-well- to double-well-potential (simulated tipping point, Fig. |3]). Our aim is to test 
if the proposed methodology is capable to capture the modelled dynamics of the series 
and forecast the record adequately. 

We perform hindcast of these series by choosing a certain point where we start 
extrapolation of the empirical probability density, then we compare the modelled series 



with actual data at the end of the forecast. Figures l|2 show datasets in panels (a), 
the probability density functions and histograms of the data and hindcast in panels 
(b,c) of each plot, and the samples of modelled data in panels (d). 

Figure |3] combines two hindcasts of the series bifurcating from one-well to double- 
well potential, in the intervals shown by arrows in the figure (probability density func- 
tions and histogram of observed and modelled data in panels (b,d,e,h) and (c,f,g,i)). 
The initial potential is U{z) = z"^ — 2z'^] then the term with 1st power of z starts 
growing gradually from zero until the potential reaches form U{z) = — 2z^ + 8z. The 
sample of the data is shown in Figure |3^. In other panels of the figure we demonstrate 
two forecasts (from 2100 and 3400). 

The forecast extends for 730 datapoint (two years of "daily" data after 800-point 
extrapolation of the probability density). In this experiment, the forecast is performed 
in total for 1500 "days" without any intermediate assimmilation of modelled data, 
and for these condition the modelled series is very close in statistical properties to 
the hindcast data (although the probability density is not entirely identical). The 
result demonstrates that for certain bifurcating systems with gradual dynamics the 
methodology may be very efiicient as a forecast tool. 

3.2. Climate data 

Similarly to artificial data in the previous section, we performed hindcast exper- 
iments with observed climate data, temperature and sea-ice extent. The results are 



shown in Figures 4|5 The datasets are shown in panels (a), the probability density 
functions and histograms of the data and hindcast in panels (b,c) of each plot, and the 
samples of modelled data in panels (d). 

The Central England Temperature [10], which is available as monthly series since 
1659 and as daily since 1772, is considered here in daily format. It is first deseasonalised 
and then the hindcast is performed as in the above artificial records. 

The fiuctuations of the Arctic sea-ice area are not only deseasonalised, but also the 
quadratic trend is removed, as we pre-processed the sea-ice data in our recent paper [9] . 
This was done to study the properties of the fiuctuations. 
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Altough the real data has much more comphcated variabihty and dynamics, the 
method performs as well as in the cases of artificial data analysed above, with modelled 
series having the same statistical properties as the initial data. 

3.3. Forecast skill 

The skill of the forecast is calculated for multiple subsets along series using the 
hindcast approach as described above. We plot the skill in Figure [6] and conclude that 
in some cases the exact statistical properties and correlations are not reproduced well, 
with acceptable mean values over subsets but rather large standard deviations. This 
is because our stochastic approach is based on the probability density function, which 
means that the exact dynamics may vary for series with the same histograms, and 
skill estimators based on point-wise comparison of time series may deviate from ideal 
values. However, the patterns of the series are very close and suitable for long-term 
stochastic predictions, as we show below for the Arctic sea-ice. 

4. Future Arctic sea-ice area dynamics 

Arctic sea-ice dynamics have been a topic of recent scientific debate, and available 
estimates of when summer ice cover will disappear range from as early as 2016 to never. 

There is also an ongoing discussion on tipping points in the Arctic region (see [131 
[HI |3] and references therein) . 

Here we propose a stochastic forecast based on the above methodology, without 
taking into account complex feedbacks (that may reverse the declining dynamics as well 
as speed it up). Assuming that the present dynamics continue its gradual development, 
and taking into account what we already know about this time series, we forecast it 
using potential analysis. 

In Fig. [7], we show a long-term forecast of the Arctic sea-ice by combining our 
knowledge about seasonality and current trends in the data. From the point of forecast 
(which is chosen in 2008 rather than most recently in order to assess the accuracy of the 
forecast over the recently observed data), we combine the last 5-year seasonal average, 
extrapolated quadratic trend, and the forecast of the fluctuations as described above. 
In the inset we show that between 2008 and 2010 the forecast is surprisingly good and 
very close to the observed data; the later departure is mainly due to very anomalous 
summers in the recent few years. 

The forecast indicates complete loss of Arctic sea-ice area by 2030s. Obviously, the 
method provides an estimate of sea ice loss in the Arctic Ocean if the system will not 
experience new feedbacks which are not included in our approach. It is therefore a 
conservative estimate solely based on extrapolation of the current trends. 
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5. Discussion. 

We have developed a forecasting technique based on the potential analysis based 
on dynamical propagation of the probability density. As the main idea, we employ ap- 
proximations of the empirical probability distribution and extrapolate them to the fu- 
ture probability distribution of data. We conducted experiments with artificial, model 
and observed climatic data and showed the efficient forecast performance of the new 
methodology. 

As one particular example, we applied the method to the expected sea ice trend in 
the Arctic Ocean. The dynamics of Arctic sea ice was forecast and shown to collaps 
in the 2030s if the current trends remain the same. This shall be compared to fully 
dynamical models. 

As a logical next step, bootstrap resampling (Section 2.4) shall help to construct 
prediction cofidence intervals based on percentiles or standard errors 
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Figure 1: Artificial data with oscillating potential, from one to two wells and back: a) time series; 
b) hindcast empirical probability density, where blue curve is the initial statistics at the beginning of 
the hindcast, black curves are extrapolated densities up to 100 time units ahead, cyan curve is the 
real pdf at the end of the forecast, for comparisotDwith extrapolation; c) histograms of the forecast 
and real data at the end of extrapolation; d) time series corresponding to histograms in the panel c). 
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Figure 2: Artificial double-well-potential data with decreasing noise level: a) time series; b) hindcast 
empirical probability density, where blue curve is the initial statistics at the beginning of the hindcast, 
black curves are extrapolated densities up to 100 time units ahead, cyan curve is the real pdf at the 
end of the forecast, for comparison with extrapolsLtton; c) histograms of the forecast and real data at 
the end of extrapolation; d) time series corresponding to histograms in the panel c). 
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Figure 3: Artificial data bifurcating from one-well to double-well dynamics, a) time series - the 
dynamics is potentially extrapolated in two intervals marked by arrows: from 2100 to 2900 and 
from 3400 to 4200; after that the series is forecast for 730 datapoints; b) potential curves for the 
interval 2100-2900: blue curve is the initial prororoility density, cyan is the final curve; red curves 
are extrapolated at equal steps; c) the same as (b) for the interval 3400-4200; (d)-(e) are histograms 
for extrapolation at point 2900: comparison of the actual histogram and the result of extrapolation; 
(f)-(g) the same as (d-e) for the point 4200; h) comparison of the forecast and actual data in interval 
2900-3630: i) comparison of the forecast and actual data in interval 4200-4930. 
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Figure 4: Central England Temperature deseasonalised fluctuations: a) time series; b) hindcast em- 
pirical probability density, where blue curve is the initial statistics at the beginning of the hindcast, 
black curves are extrapolated densities up to 100 time units ahead, cyan curve is the real pdf at the 
end of the forecast, for comparison with extrapolsLt^on; c) histograms of the forecast and real data at 
the end of extrapolation; d) time series corresponding to histograms in the panel c). 
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Figure 5: Arctic sea-ice area fluctuations after deseasonalising and removal of quadratic decreasing 
trend: a) time series; b) hindcast empirical probability density, where blue curve is the initial statistics 
at the beginning of the hindcast, black curves are extrapolated densities up to 100 time units ahead, 
cyan curve is the real pdf at the end of the foreca^ for comparison with extrapolation; c) histograms 
of the forecast and real data at the end of extrapolation; d) time series corresponding to histograms 
in the panel c). 





Figure 6: The skill statistics of the forecasts of four datasets. 
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Figure 7: Arctic sea-ice area forecast until 2035, which indicates zero level of ice in the 2030s. The 
inset plot shows magnification of the main plot for the period 2007-2015, where one can see how the 
forecast series (red) is started from 2008 and extends beyond the observed data (black). 
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